Opis konkursu¶

Zadanie polega na prognozie sprzedaży produktów z różnych dziedziń w sieci sklepów Favorita w Ekwadorze.

Zestaw treningowy i testowy zawierają następujące dane:

  • sales - tylko w zestawie treningowym, ilość sprzedanych produktów, ułamki pochodzą z produktów na wagę
  • store_nbr - numer sklepu, w sumie jest 54 różnych sklepów
  • family - dziedzina produktów np. pieczywo, nabiał, książki, etc., w sumie jest 33 różnych dziedzin
  • onpromotion - ilość produktów na promocji
  • date - data korespondująca z wyżej wymienionymi danymi danego dnia
  • id - indeks

Ponadto do wykorzystania są następujące zestawy danych:

  • stores - informacje o lokalizacji sklepów - miasto, rewir/prowincja, klasa sklepu - te informacje nie zostały wykorzystane w modelowaniu
  • oil - cena ropy danego dnia, od której gospodarka Ekwadoru jest silnie zależna
  • holiday_events - informacje dotyczące świąt i wydarzeń takie jak czy święto jest obchodzone w całym kraju, czy tylko w pewnych rejonach, czy święto/wydarzenie zostało przeniesione na inny dzień

Zestaw testowy zawiera informacje na 15 dni dla których należy przewidzieć ilości sprzedaży produktów każdej dziedziny w każdym ze sklepów.

Opis konkursu zawierał dwie dodatkowe informacje:

  • wypłaty odbywają się w połowie miesiąca
  • w połowie kwietnia 2016 roku zaszło trzęsienie ziemi) w północno zachodniej części Ekwadoru, którego efekty mogą być widoczne w sprzedaży w niektórych sklepach (zniszczona infrastruktura, pomoc charytatywna, priorytetyzowanie niezbędnych produktów)

Link do konkursu na Kaggle i oryginalnego opisu.

Spis treści:

  • Eksploracja danych
  • Analiza czasowa
  • Przygotowanie zestawów treningowych
  • Modelowanie
  • Wnioski
In [ ]:
import pandas as pd
import numpy as np

import math
import matplotlib.pyplot as plt
import seaborn as sns
import textwrap
import json

import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objs as go
import plotly.io as pio
from IPython.display import Image

from sklearn.linear_model import LinearRegression
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess

from datetime import datetime, timedelta
import warnings
In [ ]:
# Prosta funkcja do wyświetlania wykresów plotly jako obraz rastrowy:
PNG_PLOTS = True

def check_png(fig, rule=PNG_PLOTS):
    if rule:
        img_bytes = fig.to_image(format="png")
        fig = Image(img_bytes)
        return fig
    else:
        return fig
In [ ]:
df = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
holidays = pd.read_csv('holidays_events.csv')
oil = pd.read_csv('oil.csv')
stores = pd.read_csv('stores.csv')
transactions = pd.read_csv('transactions.csv')
In [ ]:
df
Out[ ]:
id date store_nbr family sales onpromotion
0 0 2013-01-01 1 AUTOMOTIVE 0.000 0
1 1 2013-01-01 1 BABY CARE 0.000 0
2 2 2013-01-01 1 BEAUTY 0.000 0
3 3 2013-01-01 1 BEVERAGES 0.000 0
4 4 2013-01-01 1 BOOKS 0.000 0
... ... ... ... ... ... ...
3000883 3000883 2017-08-15 9 POULTRY 438.133 0
3000884 3000884 2017-08-15 9 PREPARED FOODS 154.553 1
3000885 3000885 2017-08-15 9 PRODUCE 2419.729 148
3000886 3000886 2017-08-15 9 SCHOOL AND OFFICE SUPPLIES 121.000 8
3000887 3000887 2017-08-15 9 SEAFOOD 16.000 0

3000888 rows × 6 columns

Eksploracja danych

In [ ]:
# funkcja do wyświetlania informacji o ramkach danych

def df_overview(dataframe, modulo=5, head_rows=3, zip_nans=True):
  shape_rows = dataframe.shape[0] # ile wierszy
  shape_cols = dataframe.shape[1] # ile wszystkich kolumn
  shape_cols_rest =  shape_cols % modulo # ile kolumn po dzieleniu (ogon)

  main = dataframe.iloc[:,:-shape_cols_rest] # df z głównymi kolumnami
  rest = dataframe.iloc[:,-shape_cols_rest:] # df z ogonem kolumn

  main_cols = np.array(main.columns.tolist()).reshape(int(main.shape[1] / modulo), modulo) 
  # nazwy głównych kolumn w grupach
  rest_cols = rest.columns # nazwy kolumn ogon

  for i in main_cols: # do wyświetlenia kolumn z ilorazu
    head = dataframe[i].head(head_rows) # góra outputu

    uniques = [dataframe[j].nunique() for j in head.columns] # ile unikatowych wartości

    nans =  [dataframe[j].isna().sum() for j in head.columns] # ile NaN
    nans_perc = [round(dataframe[j].isna().sum() / shape_rows  * 100, 1) for j in head.columns ] # % NaN

    nans_and_perc = [(nans[i], str(nans_perc[i])+' %') for i in range(len(nans))] 
    # kombinacja obu żeby była jeden wiersz

    types = [dataframe[dataframe[j].notnull()][j].dtype for j in head.columns] # jakie typy


    if zip_nans: # czy chcemy skrócone NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans_and_perc

      head = head.rename(index={int(len(head)-3): 'Data type',
                        int(len(head)-2): 'Unique values',
                        int(len(head)-1): 'NaNs',
                        })

    else: # czy chcemy pełne NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans
      head.loc[len(head)] = nans_perc

      head = head.rename(index={int(len(head)-4): 'Data type',
                        int(len(head)-3): 'Unique values',
                        int(len(head)-2): 'NaN total',
                        int(len(head)-1): 'NaN %'
                        })

    display(head) # wyświetla kolumny z ilorazu


  for i in [rest_cols]: # do wyświetlenia kolumn z reszty, robimy [list] żeby działało tak samo
    head = dataframe[i].head(head_rows) # góra outputu

    uniques = [dataframe[j].nunique() for j in head.columns] # ile unikatowych wartości

    nans =  [dataframe[j].isna().sum() for j in head.columns] # ile NaN
    nans_perc = [round(dataframe[j].isna().sum() / shape_rows  * 100, 1) for j in head.columns ] # % NaN

    nans_and_perc = [(nans[i], str(nans_perc[i])+' %') for i in range(len(nans))] 
    # kombinacja obu żeby była jeden wiersz

    types = [dataframe[dataframe[j].notnull()][j].dtype for j in head.columns] # jakie typy
    if zip_nans: # czy chcemy skrócone NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans_and_perc
      head = head.rename(index={int(len(head)-3): 'Data type',
                        int(len(head)-2): 'Unique values',
                        int(len(head)-1): 'NaNs',
                        })

    else: # czy chcemy pełne NaN
      head.loc[len(head)] = types
      head.loc[len(head)] = uniques
      head.loc[len(head)] = nans
      head.loc[len(head)] = nans_perc
      head = head.rename(index={int(len(head)-4): 'Data type',
                        int(len(head)-3): 'Unique values',
                        int(len(head)-2): 'NaN total',
                        int(len(head)-1): 'NaN %'
                        })
    print(dataframe.shape)
    display(head) # wyświetla kolumny z reszty
In [ ]:
print('train:')
display(df_overview(df, 10))
print('test:')
display(df_overview(test, 10))
print('holidays:')
display(df_overview(holidays, 10))
print('oil:')
display(df_overview(oil, 10))
print('stores:')
display(df_overview(stores, 10))
print('transactions:')
display(df_overview(transactions, 10))
train:
(3000888, 6)
id date store_nbr family sales onpromotion
0 0 2013-01-01 1 AUTOMOTIVE 0.0 0
1 1 2013-01-01 1 BABY CARE 0.0 0
2 2 2013-01-01 1 BEAUTY 0.0 0
Data type int64 object int64 object float64 int64
Unique values 3000888 1684 54 33 379610 362
NaNs (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %)
None
test:
(28512, 5)
id date store_nbr family onpromotion
0 3000888 2017-08-16 1 AUTOMOTIVE 0
1 3000889 2017-08-16 1 BABY CARE 0
2 3000890 2017-08-16 1 BEAUTY 2
Data type int64 object int64 object int64
Unique values 28512 16 54 33 212
NaNs (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %)
None
holidays:
(350, 6)
date type locale locale_name description transferred
0 2012-03-02 Holiday Local Manta Fundacion de Manta False
1 2012-04-01 Holiday Regional Cotopaxi Provincializacion de Cotopaxi False
2 2012-04-12 Holiday Local Cuenca Fundacion de Cuenca False
Data type object object object object object bool
Unique values 312 6 3 24 103 2
NaNs (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %)
None
oil:
(1218, 2)
date dcoilwtico
0 2013-01-01 NaN
1 2013-01-02 93.14
2 2013-01-03 92.97
Data type object float64
Unique values 1218 998
NaNs (0, 0.0 %) (43, 3.5 %)
None
stores:
(54, 5)
store_nbr city state type cluster
0 1 Quito Pichincha D 13
1 2 Quito Pichincha D 13
2 3 Quito Pichincha D 8
Data type int64 object object object int64
Unique values 54 22 16 5 17
NaNs (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %) (0, 0.0 %)
None
transactions:
(83488, 3)
date store_nbr transactions
0 2013-01-01 25 770
1 2013-01-02 1 2111
2 2013-01-02 2 2358
Data type object int64 int64
Unique values 1682 54 4993
NaNs (0, 0.0 %) (0, 0.0 %) (0, 0.0 %)
None

Wykres dziedzin sprzedaży¶

In [ ]:
families_mean = df.groupby('family').mean(numeric_only=True).sort_values('sales', ascending=False)
# families_mean

fig = px.bar(x=families_mean.index, y=families_mean['sales'])
fig.update_layout(template='plotly_dark', title='Średnia dzienna sprzedaż według rodzaju', yaxis_title='Średnia sprzedaż', xaxis_title='dziedzina')
check_png(fig)
Out[ ]:
No description has been provided for this image

Wykres ilości transakcji¶

In [ ]:
fig = px.box(x=transactions['store_nbr'], y=transactions['transactions'])
fig.update_traces(jitter=0.75, marker=dict(size=3))
fig.update_layout(
    width=1500,
    template='plotly_dark', title='Wykres pudełkowy ilości transakcji w danym sklepie',
    xaxis_title="Numer sklepu",
    yaxis_title="Ilość transakcji")
check_png(fig)
Out[ ]:
No description has been provided for this image

Sklepy i dziedziny sprzedaży

In [ ]:
# Kolory z pakietu plotly
all_colors = []

for i in px.colors.qualitative.Plotly:
    all_colors.append(i)
for i in px.colors.qualitative.Set1:
    all_colors.append(i)
for i in px.colors.qualitative.Set2:
    all_colors.append(i)
for i in px.colors.qualitative.Set3:
    all_colors.append(i)
In [ ]:
stores = df.groupby(['store_nbr', 'family']).sum(numeric_only=True).reset_index()

color_map = dict(zip(stores['family'].unique(), all_colors[:33]))

# stores
# color_map
In [ ]:
fig = go.Figure()
for family, color in color_map.items():
    subset_df = stores[stores['family'] == family]
    fig.add_trace(go.Bar(
        x=subset_df['store_nbr'],
        y=subset_df['sales'],
        text=subset_df['family'],
        marker_color=color,
        name=family
    ))

fig.update_layout(
    template='plotly_dark', width=1900, height=1200,
    title='Całkowita sprzedaż w danym sklepie, zależnie od dziedziny', xaxis_title="Numer sklepu", yaxis_title="Suma sprzedaży", legend_title='Dziedzina',  
    showlegend=True,  
    barmode='stack',
    xaxis=dict(
        tickmode='linear',  
        dtick=3  
    )
)
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Promocje¶

Funckja tworząca trace dla każdej dziedziny sprzedaży i ilości promocji.

  • trace to obiekt graficzny, zawierający dane do wizualizacji wykresu
  • dane te zostają nałożone na "płótno"
  • styl wykresu można określić przy tworzeniu trace-ów lub przy nakładaniu

Do wykreślenia sprzedaży i promocji wykorzystano:

  • metodę subplots - pozwala na wyświetlenie wielu wykresów na jednej "figurze"
  • parametr secondary_y - umożliwia nałożenie na jeden wykres dwie krzywe w różnych, osobnych skalach, w tym przypadku główna (lewa) oś Y to sprzedaż produktów danej dziedziny, a drugorzędna (prawa) oś Y to ilość produktów w promocji
In [ ]:
prom_figs = []
def get_sales_and_prom(sales_data, family):
    q = sales_data.query(f'family == "{family}"')[['date', 'sales', 'onpromotion']].groupby('date').sum().reset_index()
    trace_p = go.Scatter(x=q['date'], y=q['onpromotion'], mode='markers', marker=dict(size=2))
    trace_s = go.Scatter(x=q['date'], y=q['sales'], line=dict(width=1))

    trace_set = [trace_s, trace_p, family]
    prom_figs.append(trace_set)
In [ ]:
for i in df['family'].unique():
    get_sales_and_prom(df, i)
In [ ]:
figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(figs[counter][1], row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1
print(counter)
fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)
16
Out[ ]:
No description has been provided for this image
In [ ]:
# niestety `secondary_y i szerokie `subplots` jakoś nie działało
figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 16
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(figs[counter][1], row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image
In [ ]:
fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(figs[32][0], )
fig.add_trace(figs[32][1], secondary_y=True)
fig.update_yaxes(title_text=figs[32][-1],)

fig.update_layout(height=400, width=800, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Wnioski z promocji:¶

  • w większości brak wyraźnej korelacji

  • najbardziej wyraźna w SCHOOL AND OFFICE SUPPLIES, w okolicach rozpoczęcia roku szkolnego (w Ekwadorze również początek września)

  • LAWN AND GARDEN również wykazuje duży związek między wzrostem sprzedaży i promocjami

  • podczas modelowania z lagami dla promocji wyniki były mniej trafne

Ropa / benzyna¶

  • tutaj też zostaną wykorzystane drugorzędne osie Y
  • prawdopodobnie takie rozpatrywanie nie jest konieczne, ale pozwoli na łatwe zidentyfikowanie, które dziedziny sprzedaży najbardziej korelują z ceną benzyny, tudzież ogólnym stanem gospodarki (gospodarka Ekwadoru silnie opiera się na wydobyciu ropy)
In [ ]:
oil = oil.fillna(method='backfill')
oil_c = oil[:-12]
oil_c.tail()
Out[ ]:
date dcoilwtico
1201 2017-08-09 49.59
1202 2017-08-10 48.54
1203 2017-08-11 48.81
1204 2017-08-14 47.59
1205 2017-08-15 47.57
In [ ]:
sales_per_day = df.groupby('date').sum(numeric_only=True)['sales']
In [ ]:
fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(go.Scatter(x=sales_per_day.index, y=sales_per_day, name='Total sales per day'))
fig.add_trace(go.Scatter(x=oil['date'], y=oil['dcoilwtico'], name='Daily oil price'), secondary_y=True)
fig.update_yaxes(
    
)

fig.update_layout(height=800, width=1800, template='plotly_dark',)
check_png(fig, False)
In [ ]:
figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(go.Scatter(x=oil['date'], y=oil['dcoilwtico'],), row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text='Oil', row=i, col=j, range=[-0, 120], secondary_y=True)
        counter += 1
print(counter)
fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sprzedaż i ropa/benzyna')
check_png(fig, PNG_PLOTS)
16
Out[ ]:
No description has been provided for this image
In [ ]:
figs = prom_figs
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 16
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][0], row=i, col=j, )
        fig.add_trace(go.Scatter(x=oil['date'], y=oil['dcoilwtico'], name='Daily oil price'), row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text='Oil', row=i, col=j, range=[-0, 120], secondary_y=True)
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sprzedaż i ropa/benzyna')

check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Korelacje¶

Korelacje dziedzin względem:

  • promocji
  • ropy/benzyny
In [ ]:
corrs = df.groupby(['date','family']).sum(numeric_only=True)[['sales', 'onpromotion']].unstack().corr()
corrs = corrs.reset_index()
corrs['new_family'] = corrs['level_0'] + ' ' + corrs['family']
corrs = corrs.set_index('new_family')
corrs.drop(columns=['level_0', 'family'], level=0, inplace=True)
In [ ]:
corrs = corrs.transpose().reset_index()
corrs['new_family'] = corrs['level_0'] + ' ' + corrs['family']
corrs
corrs = corrs.set_index('new_family')
corrs.drop(columns=['level_0', 'family'], inplace=True)
In [ ]:
corrs.fillna(0, inplace=True)
In [ ]:
# fig = px.imshow(corrs.round(decimals=2) * 100, # w skali 1-100 żeby zajmowało mniej miejsca w kwadracie
#                 text_auto=True,
#                 labels=dict(color="Correlation"),
#                 color_continuous_scale='RdBu_r', 
#                 color_continuous_midpoint=0)

# fig.update_traces(textfont_size=8)
# fig.update_xaxes(tickfont=dict(size=9))
# fig.update_yaxes(tickfont=dict(size=9))

# fig.update_layout(height=1200, width=1200,)

# check_png(fig)
In [ ]:
grouped = df.groupby("family")

correlations = {}

for family, group in grouped:
    correlation = group["sales"].corr(group["onpromotion"])
    correlations[family] = correlation
In [ ]:
fig = px.bar(pd.Series(correlations).sort_values(ascending=False))
fig.update_layout(template='plotly_dark', title='Korelacja sprzedaży względem promocji', 
                  width=1200, showlegend=False,
                  xaxis_title="Dziedzina sprzedaży",
                  yaxis_title="Korelacja",)
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image
In [ ]:
grouped = df.merge(right=oil_c).groupby("family")

correlations = {}

for family, group in grouped:
    correlation = group["sales"].corr(group["dcoilwtico"])
    correlations[family] = correlation
In [ ]:
fig = px.bar(pd.Series(correlations).sort_values(ascending=True))
fig.update_layout(template='plotly_dark', title='Korelacja sprzedaży względem ceny beznyny', 
                  width=1200, showlegend=False,
                  xaxis_title="Dziedzina sprzedaży",
                  yaxis_title="Korelacja",)
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Analiza czasowa

W analizie czasowej zostanie rozpatrzone jak sprzedaż zmienia się w czasie. Uwzględni dni tygodnia, średnią ruchomą, regresję liniową i wielomianową

In [ ]:
def fix_dates(df, datename='date', set_date=True):
    df['date'] = pd.to_datetime(df[datename], format='%Y-%m-%d')
    df['date'] = df.date.dt.to_period('D')
    if set_date:
        df = df.set_index('date')   
    return df
In [ ]:
grouped = df.copy()
grouped = fix_dates(grouped, 'date', False).groupby('date').mean(numeric_only=True)
In [ ]:
grouped['day of week'] = grouped.index.dayofweek

Dni tygodnia według dokumentacji pandas

In [ ]:
grouped
Out[ ]:
id store_nbr sales onpromotion day of week
date
2013-01-01 890.5 27.5 1.409438 0.000000 1
2013-01-02 2672.5 27.5 278.390807 0.000000 2
2013-01-03 4454.5 27.5 202.840197 0.000000 3
2013-01-04 6236.5 27.5 198.911154 0.000000 4
2013-01-05 8018.5 27.5 267.873244 0.000000 5
... ... ... ... ... ...
2017-08-11 2992868.5 27.5 463.733851 7.956790 4
2017-08-12 2994650.5 27.5 444.798280 4.664422 5
2017-08-13 2996432.5 27.5 485.768618 5.209315 6
2017-08-14 2998214.5 27.5 427.004717 4.513468 0
2017-08-15 2999996.5 27.5 427.980884 5.951178 1

1684 rows × 5 columns

In [ ]:
day_names = {0: 'Poniedziałek', 1: 'Wtorek', 2: 'Środa', 3: 'Czwartek', 4: 'Piątek', 5: 'Sobota', 6: 'Niedziela'}
grouped['day of week'] = grouped['day of week'].map(day_names)
In [ ]:
grouped
Out[ ]:
id store_nbr sales onpromotion day of week
date
2013-01-01 890.5 27.5 1.409438 0.000000 Wtorek
2013-01-02 2672.5 27.5 278.390807 0.000000 Środa
2013-01-03 4454.5 27.5 202.840197 0.000000 Czwartek
2013-01-04 6236.5 27.5 198.911154 0.000000 Piątek
2013-01-05 8018.5 27.5 267.873244 0.000000 Sobota
... ... ... ... ... ...
2017-08-11 2992868.5 27.5 463.733851 7.956790 Piątek
2017-08-12 2994650.5 27.5 444.798280 4.664422 Sobota
2017-08-13 2996432.5 27.5 485.768618 5.209315 Niedziela
2017-08-14 2998214.5 27.5 427.004717 4.513468 Poniedziałek
2017-08-15 2999996.5 27.5 427.980884 5.951178 Wtorek

1684 rows × 5 columns

In [ ]:
fig = px.box(x=grouped['day of week'], y=grouped['sales'])
fig.update_traces(jitter=0.75, marker=dict(size=3))
fig.update_layout(
    width=1500,
    template='plotly_dark', title='Wykres pudełkowy sprzedaży zależnie od dnia tygodnia',
    xaxis_title="Dzień tygodnia",
    yaxis_title="Sprzedaż",
    xaxis=dict(
        categoryorder="array",
        categoryarray=list(day_names.values()))
    )
check_png(fig)
Out[ ]:
No description has been provided for this image

Wnioski ze sprzedaży zależnie od dnia tygodnia:

  • największa sprzedaż w niedzielę i sobotę
  • najmniejsza w czwartek
  • zdecydowanie to rozróżnienie będzie istotne w modelowaniu

Trend¶

In [ ]:
sales_per_day = sales_per_day.to_frame()
sales_per_day['time_step'] = np.arange(len(sales_per_day.index))
sales_per_day
Out[ ]:
sales time_step
date
2013-01-01 2511.618999 0
2013-01-02 496092.417944 1
2013-01-03 361461.231124 2
2013-01-04 354459.677093 3
2013-01-05 477350.121229 4
... ... ...
2017-08-11 826373.722022 1679
2017-08-12 792630.535079 1680
2017-08-13 865639.677471 1681
2017-08-14 760922.406081 1682
2017-08-15 762661.935939 1683

1684 rows × 2 columns

In [ ]:
X = sales_per_day.loc[:, ['time_step']]
y = sales_per_day.loc[:, 'sales']
In [ ]:
model = LinearRegression()
model.fit(X, y)
Out[ ]:
LinearRegression()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
LinearRegression()
In [ ]:
y_pred_lr = pd.Series(model.predict(X), index=X.index)
In [ ]:
all_holidays = holidays.query('locale == "National"').join(sales_per_day, on='date').fillna(0)
In [ ]:
transfers = holidays.query('transferred == True')
transfers = transfers.join(sales_per_day, on='date').fillna(0)
transfers
Out[ ]:
date type locale locale_name description transferred sales time_step
19 2012-10-09 Holiday National Ecuador Independencia de Guayaquil True 0.000000 0.0
72 2013-10-09 Holiday National Ecuador Independencia de Guayaquil True 322529.418957 281.0
135 2014-10-09 Holiday National Ecuador Independencia de Guayaquil True 510519.923104 645.0
255 2016-05-24 Holiday National Ecuador Batalla de Pichincha True 606377.205216 1236.0
266 2016-07-25 Holiday Local Guayaquil Fundacion de Guayaquil True 697383.602092 1298.0
268 2016-08-10 Holiday National Ecuador Primer Grito de Independencia True 658457.436112 1314.0
297 2017-01-01 Holiday National Ecuador Primer dia del ano True 12082.500997 1457.0
303 2017-04-12 Holiday Local Cuenca Fundacion de Cuenca True 791762.312883 1558.0
312 2017-05-24 Holiday National Ecuador Batalla de Pichincha True 746303.627126 1600.0
324 2017-08-10 Holiday National Ecuador Primer Grito de Independencia True 651386.911970 1678.0
328 2017-09-28 Holiday Local Ibarra Fundacion de Ibarra True 0.000000 0.0
340 2017-12-06 Holiday Local Quito Fundacion de Quito True 0.000000 0.0
In [ ]:
window_v = 30
SMA_30 = sales_per_day['sales'].rolling(window=window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()
window_v = 365
SMA_365 = sales_per_day['sales'].rolling(window= window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()

fig = go.Figure()
fig.add_trace(go.Scatter(x=sales_per_day.index, y=sales_per_day['sales'], mode='lines', name='Total sales', 
                         line=dict(width=1.5)))
fig.add_trace(go.Scatter(x=all_holidays['date'], y=all_holidays['sales'], mode='markers', name='All holidays'))
# fig.add_trace(go.Scatter(x=transfers['date'], y=transfers['sales'], mode='markers', name='Transferred holidays'))
fig.add_trace(go.Scatter(x=y_pred_lr.index, y=y_pred_lr, mode='lines', name='Regression fit', 
                         line=dict(width=3), marker=dict(color='#FECB52')))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_30, 
                         name=f'Trend sprzedaży w przedziale czasowym 30 dni'))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_365, 
                         name=f'Trend sprzedaży w przedziale czasowym 180 dni'))


fig.update_layout(template='plotly_dark', 
                  height=700, width=1800,
                  legend=dict(
                      orientation="h",
                      yanchor="bottom",
                      y=1.02,
                      xanchor="right",
                      x=1))

check_png(fig)
Out[ ]:
No description has been provided for this image
In [ ]:
warnings.filterwarnings('ignore', category=UserWarning)
In [ ]:
average_sales = sales_per_day['sales']
y = average_sales.copy()

dp = DeterministicProcess(average_sales.index, constant=False, order=3) # do 3 potęgi, krzywa wielomianowa zamiast linii 

X = dp.in_sample()

forecast_end = 32
X_fore = dp.out_of_sample(forecast_end) # kontynuacja krzywej na 16 kroków (dni)
In [ ]:
average_sales.index
Out[ ]:
Index(['2013-01-01', '2013-01-02', '2013-01-03', '2013-01-04', '2013-01-05',
       '2013-01-06', '2013-01-07', '2013-01-08', '2013-01-09', '2013-01-10',
       ...
       '2017-08-06', '2017-08-07', '2017-08-08', '2017-08-09', '2017-08-10',
       '2017-08-11', '2017-08-12', '2017-08-13', '2017-08-14', '2017-08-15'],
      dtype='object', name='date', length=1684)
In [ ]:
X_fore_steps = X_fore.index.to_list()
print(X_fore_steps)
[1685, 1686, 1687, 1688, 1689, 1690, 1691, 1692, 1693, 1694, 1695, 1696, 1697, 1698, 1699, 1700, 1701, 1702, 1703, 1704, 1705, 1706, 1707, 1708, 1709, 1710, 1711, 1712, 1713, 1714, 1715, 1716]
In [ ]:
start_date = datetime.strptime("2017-08-16", "%Y-%m-%d")
end_date = start_date + timedelta(days = forecast_end) 
In [ ]:
datelist = pd.date_range(start="2017-08-16", end=end_date)
datelist
Out[ ]:
DatetimeIndex(['2017-08-16', '2017-08-17', '2017-08-18', '2017-08-19',
               '2017-08-20', '2017-08-21', '2017-08-22', '2017-08-23',
               '2017-08-24', '2017-08-25', '2017-08-26', '2017-08-27',
               '2017-08-28', '2017-08-29', '2017-08-30', '2017-08-31',
               '2017-09-01', '2017-09-02', '2017-09-03', '2017-09-04',
               '2017-09-05', '2017-09-06', '2017-09-07', '2017-09-08',
               '2017-09-09', '2017-09-10', '2017-09-11', '2017-09-12',
               '2017-09-13', '2017-09-14', '2017-09-15', '2017-09-16',
               '2017-09-17'],
              dtype='datetime64[ns]', freq='D')
In [ ]:
X_fore['date'] = datelist[:forecast_end]
X_fore = X_fore.set_index('date')
X_fore
Out[ ]:
trend trend_squared trend_cubed
date
2017-08-16 1685.0 2839225.0 4.784094e+09
2017-08-17 1686.0 2842596.0 4.792617e+09
2017-08-18 1687.0 2845969.0 4.801150e+09
2017-08-19 1688.0 2849344.0 4.809693e+09
2017-08-20 1689.0 2852721.0 4.818246e+09
2017-08-21 1690.0 2856100.0 4.826809e+09
2017-08-22 1691.0 2859481.0 4.835382e+09
2017-08-23 1692.0 2862864.0 4.843966e+09
2017-08-24 1693.0 2866249.0 4.852560e+09
2017-08-25 1694.0 2869636.0 4.861163e+09
2017-08-26 1695.0 2873025.0 4.869777e+09
2017-08-27 1696.0 2876416.0 4.878402e+09
2017-08-28 1697.0 2879809.0 4.887036e+09
2017-08-29 1698.0 2883204.0 4.895680e+09
2017-08-30 1699.0 2886601.0 4.904335e+09
2017-08-31 1700.0 2890000.0 4.913000e+09
2017-09-01 1701.0 2893401.0 4.921675e+09
2017-09-02 1702.0 2896804.0 4.930360e+09
2017-09-03 1703.0 2900209.0 4.939056e+09
2017-09-04 1704.0 2903616.0 4.947762e+09
2017-09-05 1705.0 2907025.0 4.956478e+09
2017-09-06 1706.0 2910436.0 4.965204e+09
2017-09-07 1707.0 2913849.0 4.973940e+09
2017-09-08 1708.0 2917264.0 4.982687e+09
2017-09-09 1709.0 2920681.0 4.991444e+09
2017-09-10 1710.0 2924100.0 5.000211e+09
2017-09-11 1711.0 2927521.0 5.008988e+09
2017-09-12 1712.0 2930944.0 5.017776e+09
2017-09-13 1713.0 2934369.0 5.026574e+09
2017-09-14 1714.0 2937796.0 5.035382e+09
2017-09-15 1715.0 2941225.0 5.044201e+09
2017-09-16 1716.0 2944656.0 5.053030e+09
In [ ]:
model = LinearRegression()
model.fit(X, y)

y_pred = pd.Series(model.predict(X), index=X.index)
y_fore = pd.Series(model.predict(X_fore), index=X_fore.index)
In [ ]:
window_v = 30
SMA_30 = sales_per_day['sales'].rolling(window=window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()
window_v = 180
SMA_365 = sales_per_day['sales'].rolling(window= window_v, 
                                     center=True, 
                                     min_periods=int(window_v/2)
                                     ).mean()

fig = go.Figure()
fig.add_trace(go.Scatter(x=sales_per_day.index, y=sales_per_day['sales'], mode='lines', name='Total sales', 
                         line=dict(width=1.5)))
fig.add_trace(go.Scatter(x=all_holidays['date'], y=all_holidays['sales'], mode='markers', name='All holidays'))
fig.add_trace(go.Scatter(x=transfers['date'], y=transfers['sales'], mode='markers', name='Transferred holidays'))
fig.add_trace(go.Scatter(x=y_pred_lr.index, y=y_pred_lr, mode='lines', name='Regression fit', 
                         line=dict(width=2), marker=dict(color='#b4913e')))
fig.add_trace(go.Scatter(x=y_pred.index, y=y_pred, mode='lines', name='Polynomial-regression fit', 
                         line=dict(width=2), marker=dict(color='#FECB52')))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_30, 
                         name=f'Trend sprzedaży w przedziale czasowym 30 dni'))
fig.add_trace(go.Scatter(x=sales_per_day.index, y=SMA_365, 
                         name=f'Trend sprzedaży w przedziale czasowym 180 dni'))
fig.add_trace(go.Scatter(x=y_fore.index, y=y_fore, 
                         name=f'Prognoza sprzedaży na {forecast_end} dni'))

fig.update_xaxes(range=[datetime.strptime("2012-12-01", "%Y-%m-%d"),
                        end_date])

fig.update_layout(template='plotly_dark', 
                  height=700, width=1800,
                  legend=dict(
                      orientation="h",
                      yanchor="bottom",
                      y=1.02,
                      xanchor="right",
                      x=1))

check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Poniżej znajdują się wykreślone zestawy innych dziedzin sprzedaży

In [ ]:
q_test = df.groupby(['date', 'family']).sum().query('family == "AUTOMOTIVE"')['sales'].to_frame()
q_test['time_step'] = np.arange(len(q_test.index))
q_test
Out[ ]:
sales time_step
date family
2013-01-01 AUTOMOTIVE 0.0 0
2013-01-02 AUTOMOTIVE 255.0 1
2013-01-03 AUTOMOTIVE 161.0 2
2013-01-04 AUTOMOTIVE 169.0 3
2013-01-05 AUTOMOTIVE 342.0 4
... ... ... ...
2017-08-11 AUTOMOTIVE 441.0 1679
2017-08-12 AUTOMOTIVE 403.0 1680
2017-08-13 AUTOMOTIVE 481.0 1681
2017-08-14 AUTOMOTIVE 292.0 1682
2017-08-15 AUTOMOTIVE 337.0 1683

1684 rows × 2 columns

In [ ]:
def get_lr(query_part):
    q_test = query_part.to_frame()
    q_test['time_step'] = np.arange(len(q_test.index))

    X = q_test.loc[:, ['time_step']]  # features
    y = q_test.loc[:, 'sales']  # target

    model = LinearRegression()
    model.fit(X, y)

    y_pred_lr = pd.Series(model.predict(X), index=X.index)
    # y_pred_lr = y_pred_lr.to_frame().reset_index().drop(columns='family').set_index('date')[0]

    return y_pred_lr

get_lr(df.groupby(['date', 'family']).sum().query('family == "AUTOMOTIVE"')['sales'])
Out[ ]:
date        family    
2013-01-01  AUTOMOTIVE    241.322651
2013-01-02  AUTOMOTIVE    241.427398
2013-01-03  AUTOMOTIVE    241.532144
2013-01-04  AUTOMOTIVE    241.636891
2013-01-05  AUTOMOTIVE    241.741637
                             ...    
2017-08-11  AUTOMOTIVE    417.191855
2017-08-12  AUTOMOTIVE    417.296601
2017-08-13  AUTOMOTIVE    417.401348
2017-08-14  AUTOMOTIVE    417.506094
2017-08-15  AUTOMOTIVE    417.610840
Length: 1684, dtype: float64
In [ ]:
def get_SMA(window, sales):
    SMA = sales.rolling(window=window, 
                                     center=True, 
                                     min_periods=int(window/2)
                                     ).mean()
    return SMA

get_SMA(30, df.query(f"family == 'AUTOMOTIVE'").groupby('date').sum(numeric_only=True)['sales'])
Out[ ]:
date
2013-01-01    208.533333
2013-01-02    206.187500
2013-01-03    205.352941
2013-01-04    203.388889
2013-01-05    207.473684
                 ...    
2017-08-11    406.700000
2017-08-12    410.263158
2017-08-13    413.944444
2017-08-14    406.470588
2017-08-15    398.812500
Name: sales, Length: 1684, dtype: float64
In [ ]:
def get_dt_forecast(start_date=start_date, end_date=end_date, datelist=datelist, sales=None):
    y = sales.copy()  # the target
    dp = DeterministicProcess(average_sales.index, constant=False, order=3)
    X = dp.in_sample()
    X_fore = dp.out_of_sample(forecast_end)
    X_fore_steps = X_fore.index.to_list()

    X_fore['date'] = datelist[:forecast_end]
    X_fore = X_fore.set_index('date')

    model = LinearRegression()
    model.fit(X, y)

    y_pred = pd.Series(model.predict(X), index=X.index)
    y_fore = pd.Series(model.predict(X_fore), index=X_fore.index)

    return [y_pred, y_fore]

# get_dt_forecast(30, sales = df.query(f"family == 'AUTOMOTIVE'").groupby('date').sum(numeric_only=True)['sales']) 
In [ ]:
figs = []
for i in df['family'].unique():
    q = df.query(f"family == '{i}'").groupby('date').sum(numeric_only=True)['sales']
    trace_base = go.Scatter(x=q.index, y=q, name=f'Sprzedaż {i}', line=dict(width=1), marker=dict(color=px.colors.qualitative.Plotly[0]))

    lr = get_lr(q)
    trace_lr = go.Scatter(x=lr.index, y=lr, name=f'Regresja liniowa {i}', line=dict(width=2), marker=dict(color='#b4913e'))

    dt_forecast = get_dt_forecast(sales=q)
    trace_dt = go.Scatter(x=dt_forecast[0].index, y=dt_forecast[0], name=f'Trend {i}', line=dict(width=1.5), marker=dict(color='#FECB52'))
    trace_fr = go.Scatter(x=dt_forecast[1].index, y=dt_forecast[1], name=f'Prognoza {i}', line=dict(width=1.5), marker=dict(color=px.colors.qualitative.Plotly[2]))

    c_SMA_30 = get_SMA(30, q)
    c_SMA_180 = get_SMA(180, q)
    trace_SMA_30 = go.Scatter(x=c_SMA_30.index, y=c_SMA_30, name=f'Trend sprzedaży w przedziale czasowym 30 dni', marker=dict(color=px.colors.qualitative.Plotly[6]))
    trace_SMA_180 = go.Scatter(x=c_SMA_180.index, y=c_SMA_180, name=f'Trend sprzedaży w przedziale czasowym 180 dni', marker=dict(color=px.colors.qualitative.Plotly[4]))



    trace_set = []
    trace_set.append(trace_base)
    trace_set.append(trace_lr)
    trace_set.append(trace_dt)
    trace_set.append(trace_fr)
    trace_set.append(trace_SMA_30)
    trace_set.append(trace_SMA_180)

    trace_set.append(i)
    figs.append(trace_set)
In [ ]:
fig = make_subplots(rows=4, cols=4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        for k in range(0,6):
            fig.add_trace(figs[counter][k], row=i, col=j, )
            fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1


fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)

                  

check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image
In [ ]:
fig = make_subplots(
  rows=5, cols=4,
  specs=[
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {"rowspan":1, "colspan": 2}, None,],
      [{"rowspan":1, "colspan": 2}, None, {"rowspan":1, "colspan": 2}, None]
         ],
  print_grid=False)

counter = 16
for i in range(1,4):
    for j in range(1,5):
        for k in range(0,6):
            fig.add_trace(figs[counter][k], row=i, col=j, )
            fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        counter += 1

for j in range(1,4):
    for k in range(0,6):
        fig.add_trace(figs[counter][k], row=4, col=j, )
        fig.update_yaxes(title_text=figs[counter][-1], row=4, col=j, )
    counter += 1

for k in range(0,6):
        fig.add_trace(figs[counter][k], row=5, col=1, )
        fig.update_yaxes(title_text=figs[counter][-1], row=5, col=1, )
counter += 1

for k in range(0,6):
        fig.add_trace(figs[counter][k], row=5, col=3, )
        fig.update_yaxes(title_text=figs[counter][-1], row=5, col=3, )

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False)
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Wnioski z trendu:¶

  • większość sprzedaży stopniowo rośnie z czasem
  • koniec 2017 roku często wskazuje na spadek
  • w wielu przypadkach mamy zupełne braki sprzedaży, które wprowadzają błędy do analizy - modelowanie zostanie przeprowadzona na podzbiorach przefiltrowanych do opdowiednich dat

Dziedziny bez zarzutu:

['AUTOMOTIVE', 'BEAUTY', 'BREAD/BAKERY', 'CLEANING', 'DELI', 'EGGS', 'GROCERY I', 'GROCERY II', 'HARDWARE',]

['HOME APPLIANCES', 'LINGERIE', 'PERSONAL CARE', 'PREPARED FOODS', 'SEAFOOD']

Dziedziny dobre, ale do ucięcia w czasie:

  • BABY CARE - od 2014-03-1 i edytować wartość z 2015-08-4 lub lepiej zacząć od 2015-12-7, sprawdzić onpromotion
  • BEVERAGES - od 2015-05-28
  • CELEBRATION - od 2015-06-01
  • DIARY - od 2013-11-7
  • FROZEN FOODS - może normalnie, sprawdzić dokąd sięgają dane zestawu testowego, może zrobić one-hot encoding dla grudnia?
  • HOME AND KITCHEN I - od 2014-09-1
  • HOME AND KITCHEN II - od 2015-10-31
  • HOME CARE - od 2015-05-4
  • LADIESWEAR - od 2015-06-1
  • LAWN AND GARDEN - od 2016-12-3, może zmienić wartości dla 2017-02-13 -14 i 2017-05-12 -14
  • LIQUOR,WINE,BEER - od 2016-05-10 - wcześniej alkohol nie mógł byś sprzedawany w niedziele, od mniej więcej tego momentu już może być sprzedawany w każdy dzień tygodnia
  • MAGAZINES - od 2015-09-30
  • MEATS - edytować dla 2016-10-7 lub sprawdzić onpromotion
  • PET SUPPLIES - od 2015-06-1
  • PLAYERS AND ELECTRONICS - od 2015-06-1
  • POULTRY - od 2013-11-4
  • PRODUCE - od 2015-06-1

Słabe / niepewne dziedziny:

  • BOOKS - może od 2017-01-1
  • SCHOOLS AND OFFICE SUPPLIES - może wprowadzić wydarzenie początek roku szkolnego i może lag dla niego?

Trend - ropa/benzyna¶

In [ ]:
oil.head()
Out[ ]:
date dcoilwtico
0 2013-01-01 93.14
1 2013-01-02 93.14
2 2013-01-03 92.97
3 2013-01-04 93.12
4 2013-01-07 93.20
In [ ]:
oil_SMA30 = get_SMA(30, oil.set_index('date')['dcoilwtico'])
oil_SMA30_trace = go.Scatter(x=oil_SMA30.index, y=oil_SMA30, name ='Trend ceny ropy')
In [ ]:
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][4], row=i, col=j, )
        fig.add_trace(oil_SMA30_trace, row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text="OIL SMA", row=i, col=j, secondary_y=True)
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Średnia ruchoma dla 30 dni')

check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image
In [ ]:
fig = make_subplots(rows=4, cols=4, specs=[[{"secondary_y": True}]*4]*4)
counter = 16
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(figs[counter][4], row=i, col=j, )
        fig.add_trace(oil_SMA30_trace, row=i, col=j, secondary_y=True)
        fig.update_yaxes(title_text=figs[counter][-1], row=i, col=j, )
        fig.update_yaxes(title_text="OIL SMA", row=i, col=j, secondary_y=True)
        counter += 1

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Średnia ruchoma dla 30 dni',)

check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image
In [ ]:
fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(figs[32][4],)
fig.add_trace(oil_SMA30_trace, secondary_y=True, )
fig.update_yaxes(title_text=figs[32][-1],)
fig.update_yaxes(title_text="OIL SMA", secondary_y=True)

fig.update_layout(template='plotly_dark', 
                  height=400, width=1200,
                  legend=dict(
                      orientation="h",
                      yanchor="bottom",
                      y=1.02,
                      xanchor="right",
                      x=1))
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Sezonowość¶

  • dzień tygodnia zostanie zakodowany one-hot (indicators)
  • sezonowość zostanie zbadana poprzez periodogram
In [ ]:
def trace_periodogram(ts, detrend='linear', name='Periodogram'):
    from scipy.signal import periodogram
    fs = pd.Timedelta("365D") / pd.Timedelta("1D")
    freqencies, spectrum = periodogram(
        ts,
        fs=fs,
        detrend=detrend,
        window="boxcar",
        scaling='spectrum',
    )
    df = pd.DataFrame({'Frequency': freqencies, 'Spectrum': spectrum})
    
    trace = go.Scatter(
        x=df['Frequency'],
        y=df['Spectrum'],
        mode='lines',
        name=name,
    )
    return trace
In [ ]:
periodograms = []
for i in df['family'].unique():
    q = df.query(f'family == "{i}"').groupby('date').sum(numeric_only=True)['sales']
    periodograms.append([trace_periodogram(q, name=i), i])

periodograms = [[trace, 'Var: ' + title] for trace, title in periodograms]
In [ ]:
fig = make_subplots(rows=4, cols=4)
counter_a = 0
for i in range(1,5):
    for j in range(1,5):
        fig.add_trace(periodograms[counter_a][0], row=i, col=j, )
        fig.update_yaxes(title_text=periodograms[counter_a][-1], row=i, col=j, )  
        counter_a += 1

fig.update_xaxes(tickangle=30, tickvals=[1, 2, 4, 6, 12, 26, 52, 104],
        ticktext=[
            "Annual (1)",
            "Semiannual (2)",
            "Quarterly (4)",
            "Bimonthly (6)",
            "Monthly (12)",
            "Biweekly (26)",
            "Weekly (52)",
            "Semiweekly (104)",],)

fig.update_xaxes(type="log")
fig.update_layout(width=1500)
fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sezonowość - wariancja sprzedaży w zależności od przedziału czasowego')
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image
In [ ]:
fig = make_subplots(
  rows=5, cols=4,
  specs=[
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {}, {}],
      [{}, {}, {"rowspan":1, "colspan": 2}, None,],
      [{"rowspan":1, "colspan": 2}, None, {"rowspan":1, "colspan": 2}, None]
         ],

  print_grid=False)

counter_b = 16

for i in range(1,4):
    for j in range(1,5):
        fig.add_trace(periodograms[counter_b][0], row=i, col=j, )
        fig.update_yaxes(title_text=periodograms[counter_b][-1], row=i, col=j, )  
        counter_b += 1

for j in range(1,4):
    fig.add_trace(periodograms[counter_b][0], row=4, col=j, )
    fig.update_yaxes(title_text=periodograms[counter_b][-1], row=4, col=j, )  
    counter_b += 1

fig.add_trace(periodograms[counter_b][0], row=5, col=1, )
fig.update_yaxes(title_text=periodograms[counter_b][-1], row=5, col=1, )  
counter_b += 1
fig.add_trace(periodograms[counter_b][0], row=5, col=3, )
fig.update_yaxes(title_text=periodograms[counter_b][-1], row=5, col=3, )  

fig.update_xaxes(
        tickvals=[1, 2, 4, 6, 12, 26, 52, 104],
        ticktext=[
            "Rocznie (1)",
            "Półrocznie (2)",
            "Kwartalnie (4)",
            "Co 2 miesiące (6)",
            "Miesięcznie (12)",
            "Co 2 tygodnie (26)",
            "Tygodniowo (52)",
            "x2 w tygodniu (104)",
        ],
        tickangle=30,
    )
fig.update_xaxes(type="log")
fig.update_layout(width=1500)

fig.update_layout(height=1200, width=2200, template='plotly_dark', showlegend=False,
                  title='Sezonowość - wariancja sprzedaży w zależności od przedziału czasowego')
check_png(fig, PNG_PLOTS)
Out[ ]:
No description has been provided for this image

Sezonowość wnioski¶

(wnioski również oparte na ogólnej historii sprzedaży)

  • wiekszość najczęściej wybieranych produktów, czyli najczęściej artykuły spożywcze - GROCERY I, MEATS, SEAFOOD etc. mają bardzo dobrze widoczną sezonowość tygodniową

    • duża wariancja wewnątrz tygodnia - najmniejsze zakupy w czwartki, największe w niedziele jak już wcześniej stwierdzono
    • tu będzie one-hot encoding dni tygodnia czyli tzw. time indicators
  • dziedziny w których widać największą sezonowość w dużych okresach czasu to:

    • FROZEN FOOD - ogromny skok przed Świętami Bożego Narodzenia
    • SCHOOL AND OFFICE SUPPLIES - tu głównie wrzesień - rok szkolny zaczyna się we wrześniu, ale mamy też wzrost sprzedaży w maju/kwietniu 2016 i 2017 roku
  • duża sezonowość w skali ponad roku wynika też z okresów gdzie sprzedaż nie odbywała się, lub sprzedaż tych produktów nie była rejestrowana albo była inaczej kwalifikowana

    • przykładowo BABY CARE ma niemal zerową sprzedaż do początku 2014 roku, oraz zupełny brak sprzedaż od stycznia do maja 2015 roku
    • podobnie PET SUPPLIES
  • przykładem innego przypadku może być drób: POULTRY - do końca października 2014 roku ilość sprzedaży była około dwa razy mniejsza niż po październiku 2014 roku, przypuszczam że taki gwałtowny wzrost mógłby być spowodowany przez:

    • rozpoczęcie sprzedaży w nowych sklepach / zmianę typu sklepu
    • zmianę w cenie produktów, które np. dotychczas miały zaniżoną cenę lub zaszła zmiana w metodzie ich opodatkowania - drób i produkty mleczne są "codziennymi" typami produktów, bez których większość ludzi nie może się obejść (w szczególności drób jest w większości tańszym mięsem)
    • zmianą strategii firmy - np. dana firma, z której pochodzi ten zestaw danych, zaczęła sprowadzać więcej tych produktów do Ekwadoru, przez co ich sprzedaż wzrosła

Przygotowanie zestawów treningowych

  • ta cześć notesu wykonywana w osobno od eksploracji danych
  • idea przygotowania polega na rozdzieleniu pierwotnego zestawu testowego na 54 podzbiory dla każdego ze sklepów
  • przed podziałem utworzono dodatkowe cechy, takie jak dzień tygodnia, miesiąc, rok, dzień miesiąca itd. oraz trend i cykle Fouriera
  • oil oraz święta będą dołączane przy każdym podziale podzbioru na dziedziny i prognozowaniu

Innymi słowy, prognozowanie będzie odbywało się dla każdej osobnej dziedziny w każdym osobnym sklepie, czyli zostanie wykonane $33 \times 54 = 1782$ osobnych prognoz

In [ ]:
# czy przeprowadzić zapis w tym pokazowym zeszycie
zapis = False

Holiday dummies¶

  • dummies świąt to zero-jedynkowo zapisana tabela, w której każde święto stanowi osobną kolumnę, a wiersze reprezentowane są przez daty
  • święta i wydarzenia zostały odfiltrowane do tych, które nie zostały przeniesione i tych które są obchodzone w całym kraju
In [ ]:
train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
oil = pd.read_csv('oil.csv')
holidays = pd.read_csv('holidays_events.csv')
In [ ]:
train = fix_dates(train)
test = fix_dates(test)
oil = fix_dates(oil)
holidays = fix_dates(holidays)
In [ ]:
oil
Out[ ]:
dcoilwtico
date
2013-01-01 NaN
2013-01-02 93.14
2013-01-03 92.97
2013-01-04 93.12
2013-01-07 93.20
... ...
2017-08-25 47.65
2017-08-28 46.40
2017-08-29 46.46
2017-08-30 45.96
2017-08-31 47.26

1218 rows × 1 columns

In [ ]:
holidays
Out[ ]:
type locale locale_name description transferred
date
2012-03-02 Holiday Local Manta Fundacion de Manta False
2012-04-01 Holiday Regional Cotopaxi Provincializacion de Cotopaxi False
2012-04-12 Holiday Local Cuenca Fundacion de Cuenca False
2012-04-14 Holiday Local Libertad Cantonizacion de Libertad False
2012-04-21 Holiday Local Riobamba Cantonizacion de Riobamba False
... ... ... ... ... ...
2017-12-22 Additional National Ecuador Navidad-3 False
2017-12-23 Additional National Ecuador Navidad-2 False
2017-12-24 Additional National Ecuador Navidad-1 False
2017-12-25 Holiday National Ecuador Navidad False
2017-12-26 Additional National Ecuador Navidad+1 False

350 rows × 5 columns

In [ ]:
dum_holidays = pd.get_dummies(holidays.query('transferred == @False & locale == "National"')[['description']], drop_first=False).astype('uint8')
In [ ]:
dum_holidays.head(3)
Out[ ]:
description_Batalla de Pichincha description_Black Friday description_Carnaval description_Cyber Monday description_Dia de Difuntos description_Dia de la Madre description_Dia de la Madre-1 description_Dia del Trabajo description_Inauguracion Mundial de futbol Brasil description_Independencia de Cuenca ... description_Terremoto Manabi+5 description_Terremoto Manabi+6 description_Terremoto Manabi+7 description_Terremoto Manabi+8 description_Terremoto Manabi+9 description_Traslado Batalla de Pichincha description_Traslado Independencia de Guayaquil description_Traslado Primer Grito de Independencia description_Traslado Primer dia del ano description_Viernes Santo
date
2012-08-10 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
2012-10-12 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 1 0 0 0
2012-11-02 0 0 0 0 1 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0

3 rows × 72 columns

In [ ]:
if zapis:
    dum_holidays.to_csv('holiday_dummies.csv', index=True)

Feature engineering¶

  • tu zostaną dodane wcześniej wymienione cechy
Dzień tygodnia¶
In [ ]:
train['day_of_week'] = train.index.dayofweek
day_names = {0: 'Monday', 1: 'Tuesday', 2: 'Wednesday', 3: 'Thursday', 4: 'Friday', 5: 'Saturday', 6: 'Sunday'}
train['day_of_week'] = train['day_of_week'].map(day_names)

dum_day = pd.get_dummies(train["day_of_week"], drop_first=False)
train = pd.concat([train, dum_day], axis=1)
train.drop("day_of_week", axis=1, inplace=True)
In [ ]:
test['day_of_week'] = test.index.dayofweek
day_names = {0: 'Monday', 1: 'Tuesday', 2: 'Wednesday', 3: 'Thursday', 4: 'Friday', 5: 'Saturday', 6: 'Sunday'}
test['day_of_week'] = test['day_of_week'].map(day_names)

dum_day = pd.get_dummies(test["day_of_week"], drop_first=False)
test = pd.concat([test, dum_day], axis=1)
test.drop("day_of_week", axis=1, inplace=True)
Inne okresy¶
  • kwartał, miesiąc, rok, dzień w roku, etc.
In [ ]:
def get_date_details(df):
    df['quarter'] = df.index.quarter
    df['month'] = df.index.month
    df['year'] = df.index.year
    df['dayofyear'] = df.index.dayofyear
    df['dayofmonth'] = df.index.day
    df['weekofyear'] = df.index.week
    return df
In [ ]:
train = get_date_details(train)
test = get_date_details(test)
In [ ]:
test.head()
Out[ ]:
id store_nbr family onpromotion Friday Monday Saturday Sunday Thursday Tuesday Wednesday quarter month year dayofyear dayofmonth weekofyear
date
2017-08-16 3000888 1 AUTOMOTIVE 0 0 0 0 0 0 0 1 3 8 2017 228 16 33
2017-08-16 3000889 1 BABY CARE 0 0 0 0 0 0 0 1 3 8 2017 228 16 33
2017-08-16 3000890 1 BEAUTY 2 0 0 0 0 0 0 1 3 8 2017 228 16 33
2017-08-16 3000891 1 BEVERAGES 20 0 0 0 0 0 0 1 3 8 2017 228 16 33
2017-08-16 3000892 1 BOOKS 0 0 0 0 0 0 0 1 3 8 2017 228 16 33
Ropa/benzyna¶
In [ ]:
oil.fillna(method='bfill', inplace=True) # w niektóre dni brakuje informacji o cenie - są NaN
In [ ]:
train = train.join(oil.fillna(method='bfill'), how='left').fillna(method='bfill') # dołączenie tabelki do zestawów, uzupełnienie brakujących danych - nie ma niektórych dni w tabeli `oil`
test = test.join(oil.fillna(method='bfill'), how='left').fillna(method='bfill')

Cykle Fouriera¶

  • cykle fouriera pozwalają na wielokrotnie oszczędniejszy zapis cykliczności niż one-hot encoding każdej daty / kroku
  • opiera się na parach krzywych sinusoid i cosinusoid o określonej "długości", odwzorowywujące sezonowość danych wartości
  • ta część również była wykonywana w osobnym notesie

Optymalne typy danych przed zapisem podzbiorów

In [ ]:
train['id'] = train['id'].astype('uint32')
train['store_nbr'] = train['store_nbr'].astype('uint8')
train['family'] = train['family'].astype('category')
train['onpromotion'] = train['onpromotion'].astype('uint16')

train['quarter'] = train['quarter'].astype('uint8')
train['month'] = train['month'].astype('uint8')
train['year'] = train['year'].astype('uint16')
train['dayofyear'] = train['dayofyear'].astype('uint16')
train['dayofmonth'] = train['dayofmonth'].astype('uint8')
train['weekofyear'] = train['weekofyear'].astype('uint8')

test['id'] = test['id'].astype('uint32')
test['store_nbr'] = test['store_nbr'].astype('uint8')
test['family'] = test['family'].astype('category')
test['onpromotion'] = test['onpromotion'].astype('uint16')

test['quarter'] = test['quarter'].astype('uint8')
test['month'] = test['month'].astype('uint8')
test['year'] = test['year'].astype('uint16')
test['dayofyear'] = test['dayofyear'].astype('uint16')
test['dayofmonth'] = test['dayofmonth'].astype('uint8')
test['weekofyear'] = test['weekofyear'].astype('uint8')
In [ ]:
fourier = CalendarFourier(freq="A", order=10)  # 10 par sin/cos dla corocznej sezonowości ("A" - Annual)
dp = DeterministicProcess(
    index=train.index,
    constant=True,               # cecha dla wyrazu wolnego
    order=1,                     # trend (order=1 oznacza regresję liniową)
    seasonal=True,               # sezonowość tygodniowa (indicators)
    additional_terms=[fourier],  # sezonowość w ciągu roku (fourier)
    drop=True,                   
)

fourier_in_sample = dp.in_sample()  
fourier_in_sample
Out[ ]:
const trend s(2,7) s(3,7) s(4,7) s(5,7) s(6,7) s(7,7) sin(1,freq=A-DEC) cos(1,freq=A-DEC) ... sin(6,freq=A-DEC) cos(6,freq=A-DEC) sin(7,freq=A-DEC) cos(7,freq=A-DEC) sin(8,freq=A-DEC) cos(8,freq=A-DEC) sin(9,freq=A-DEC) cos(9,freq=A-DEC) sin(10,freq=A-DEC) cos(10,freq=A-DEC)
date
2013-01-01 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 1.000000 ... 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000
2013-01-01 1.0 2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.000000 1.000000 ... 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000
2013-01-01 1.0 3.0 0.0 1.0 0.0 0.0 0.0 0.0 0.000000 1.000000 ... 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000
2013-01-01 1.0 4.0 0.0 0.0 1.0 0.0 0.0 0.0 0.000000 1.000000 ... 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000
2013-01-01 1.0 5.0 0.0 0.0 0.0 1.0 0.0 0.0 0.000000 1.000000 ... 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2017-08-15 1.0 3000884.0 0.0 0.0 0.0 1.0 0.0 0.0 -0.680773 -0.732494 ... -0.976011 -0.217723 0.863142 -0.504961 -0.288482 0.957485 -0.440519 -0.897743 0.933837 0.357698
2017-08-15 1.0 3000885.0 0.0 0.0 0.0 0.0 1.0 0.0 -0.680773 -0.732494 ... -0.976011 -0.217723 0.863142 -0.504961 -0.288482 0.957485 -0.440519 -0.897743 0.933837 0.357698
2017-08-15 1.0 3000886.0 0.0 0.0 0.0 0.0 0.0 1.0 -0.680773 -0.732494 ... -0.976011 -0.217723 0.863142 -0.504961 -0.288482 0.957485 -0.440519 -0.897743 0.933837 0.357698
2017-08-15 1.0 3000887.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.680773 -0.732494 ... -0.976011 -0.217723 0.863142 -0.504961 -0.288482 0.957485 -0.440519 -0.897743 0.933837 0.357698
2017-08-15 1.0 3000888.0 1.0 0.0 0.0 0.0 0.0 0.0 -0.680773 -0.732494 ... -0.976011 -0.217723 0.863142 -0.504961 -0.288482 0.957485 -0.440519 -0.897743 0.933837 0.357698

3000888 rows × 28 columns

Podział na podzbiory¶

In [ ]:
train = pd.read_csv('train w features.csv',
                 dtype = {'id': 'uint32', 
                          'store_nbr': 'uint8', 
                          'family': 'category',
                          'onpromotion': 'uint16',

                          'quarter': 'uint8', 
                          'month': 'uint8', 
                          'year': 'uint16',
                          'dayofyear': 'uint16',
                          'dayofmonth': 'uint8',
                          'weekofyear': 'uint8',

                          'Friday' : 'uint8', 
                          'Monday' : 'uint8', 
                          'Saturday' : 'uint8', 
                          'Sunday' : 'uint8', 
                          'Thursday' : 'uint8', 
                          'Tuesday' : 'uint8', 
                          'Wednesday' : 'uint8', 
                          
                          'dcoilwtico' : 'float64',

                          'const': 'uint8', 
                          'trend': 'uint16',
                          })

test = pd.read_csv('test w features.csv',
                 dtype = {'id': 'uint32', 
                          'store_nbr': 'uint8', 
                          'family': 'category',
                          'onpromotion': 'uint16',

                          'quarter': 'uint8', 
                          'month': 'uint8', 
                          'year': 'uint16',
                          'dayofyear': 'uint16',
                          'dayofmonth': 'uint8',
                          'weekofyear': 'uint8',

                          'Friday' : 'uint8', 
                          'Monday' : 'uint8', 
                          'Saturday' : 'uint8', 
                          'Sunday' : 'uint8', 
                          'Thursday' : 'uint8', 
                          'Tuesday' : 'uint8', 
                          'Wednesday' : 'uint8', 
                          
                          'dcoilwtico' : 'float64',

                          'const': 'uint8', 
                          'trend': 'uint16',
                          })

# oil = pd.read_csv('oil ready.csv')
# holiday_dummies = pd.read_csv('holiday_dummies.csv').fillna(0)
In [ ]:
print(train.shape, test.shape)
(3000888, 42) (28512, 41)
In [ ]:
test = test.rename(columns={'Unnamed: 0' : 'date'})
In [ ]:
print(train.dtypes)
# print(test.dtypes)
date                    object
id                      uint32
store_nbr                uint8
family                category
sales                  float64
onpromotion             uint16
Friday                   uint8
Monday                   uint8
Saturday                 uint8
Sunday                   uint8
Thursday                 uint8
Tuesday                  uint8
Wednesday                uint8
quarter                  uint8
month                    uint8
year                    uint16
dayofyear               uint16
dayofmonth               uint8
weekofyear               uint8
dcoilwtico             float64
const                    uint8
trend                   uint16
sin(1,freq=A-DEC)      float64
cos(1,freq=A-DEC)      float64
sin(2,freq=A-DEC)      float64
cos(2,freq=A-DEC)      float64
sin(3,freq=A-DEC)      float64
cos(3,freq=A-DEC)      float64
sin(4,freq=A-DEC)      float64
cos(4,freq=A-DEC)      float64
sin(5,freq=A-DEC)      float64
cos(5,freq=A-DEC)      float64
sin(6,freq=A-DEC)      float64
cos(6,freq=A-DEC)      float64
sin(7,freq=A-DEC)      float64
cos(7,freq=A-DEC)      float64
sin(8,freq=A-DEC)      float64
cos(8,freq=A-DEC)      float64
sin(9,freq=A-DEC)      float64
cos(9,freq=A-DEC)      float64
sin(10,freq=A-DEC)     float64
cos(10,freq=A-DEC)     float64
dtype: object
In [ ]:
def fix_dates(df, datename, set_date=True):
    df['date'] = pd.to_datetime(df[datename], format='%Y-%m-%d')
    df['date'] = df.date.dt.to_period('D')
    if set_date:
        df = df.set_index('date')   
    return df

train = fix_dates(train, 'date')
test = fix_dates(test, 'date')
In [ ]:
test.head()
Out[ ]:
id store_nbr family onpromotion Friday Monday Saturday Sunday Thursday Tuesday ... sin(6,freq=A-DEC) cos(6,freq=A-DEC) sin(7,freq=A-DEC) cos(7,freq=A-DEC) sin(8,freq=A-DEC) cos(8,freq=A-DEC) sin(9,freq=A-DEC) cos(9,freq=A-DEC) sin(10,freq=A-DEC) cos(10,freq=A-DEC)
date
2017-08-16 3000888 1 AUTOMOTIVE 0 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452
2017-08-16 3000889 1 BABY CARE 0 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452
2017-08-16 3000890 1 BEAUTY 2 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452
2017-08-16 3000891 1 BEVERAGES 20 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452
2017-08-16 3000892 1 BOOKS 0 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452

5 rows × 40 columns

In [ ]:
families = test['family'].unique()
stores = test['store_nbr'].unique()
In [ ]:
for i in families[:5]:
    print(i)
AUTOMOTIVE
BABY CARE
BEAUTY
BEVERAGES
BOOKS
In [ ]:
display(
    train.query('family == @families[1] & store_nbr == @stores[0]').head()
)
id store_nbr family sales onpromotion Friday Monday Saturday Sunday Thursday ... sin(6,freq=A-DEC) cos(6,freq=A-DEC) sin(7,freq=A-DEC) cos(7,freq=A-DEC) sin(8,freq=A-DEC) cos(8,freq=A-DEC) sin(9,freq=A-DEC) cos(9,freq=A-DEC) sin(10,freq=A-DEC) cos(10,freq=A-DEC)
date
2013-01-01 1 1 BABY CARE 0.0 0 0 0 0 0 0 ... 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000 0.000000 1.000000
2013-01-02 1783 1 BABY CARE 0.0 0 0 0 0 0 0 ... 0.103102 0.994671 0.120208 0.992749 0.137279 0.990532 0.154309 0.988023 0.171293 0.985220
2013-01-03 3565 1 BABY CARE 0.0 0 0 0 0 0 1 ... 0.205104 0.978740 0.238673 0.971100 0.271958 0.962309 0.304921 0.952378 0.337523 0.941317
2013-01-04 5347 1 BABY CARE 0.0 0 1 0 0 0 0 ... 0.304921 0.952378 0.353676 0.935368 0.401488 0.915864 0.448229 0.893919 0.493776 0.869589
2013-01-05 7129 1 BABY CARE 0.0 0 0 0 1 0 0 ... 0.401488 0.915864 0.463550 0.886071 0.523416 0.852078 0.580800 0.814046 0.635432 0.772157

5 rows × 41 columns

In [ ]:
if zapis:
    for i in train['store_nbr'].unique():
        print(i,'-shape: ', {train.query('store_nbr == @i').shape})
        train.query('store_nbr == @i').to_csv(f'subsets/store_nbr_{i}.csv', index=True) # zapis podzbiorów

Modelowanie

  • prognoza będzie polegała na kombinacji dwóch modeli:
    • prostym np. regresja liniowa oparta na celu (sales) i kroku (1 - pierwszy dzień, 1688 - ostatni dzień w przypadku gdy brane byłyby wartości od początku zestawu testowego czyli od 1-01-2013)
    • złożonym np. niezawodny XGBoost oparty na wszystkich pozostałych cechach zestawu
  • od faktycznych wartości sprzedaży zostaną usunięte wartości z prostej z regresji liniowej co w wyniku da zdetrendowane wartości sprzedaży, na których będzie uczył się i przewidywał model złożony

Prognozowanie odbywało się osobno dla każdej dziedziny w każdym sklepie, czyli niecałe 1800 modeli, zajmowało to od 12 do 15 minut w zależności od hiperparametrów i złożoności modeli

Sprawdzono kilka sposobów filtrowania dat, ostatecznie zastosowano dwie metody jednocześnie:

  • filtr na podstawie eksploracji danych np. ograniczenie zestawu treningowego do wartości tylko po brakach 2015 roku
  • filtr na podstawie sumy kumulacyjnej:
    • w przypadkach gdy pierwszy filtr pozostawił daty bez sprzedaży (np. pierwszy filtr mamy od 2015-06-01), a w podzbiorze sprzedaż zaczyna się od 2016-01-01 i mamy kilka miesięcy braków
    • przy pomocy sumy kumulacyjnej można sprawdzić przy której dacie zaczyna się rzeczywista sprzedaż - dla powyższego przykładu od 2015-06-01 do 2015-12-31 suma kumulacyjna wynosiłaby 0, a w 2016-01-01 już np. 40
    • kolejno przy pomocy maski wartości True/False czy wartość jest nierówna 0 można odfiltrować puste daty

Wczytanie danych z optymalnym formatem:

In [ ]:
test = pd.read_csv('test w features.csv',
                 dtype = {'id': 'uint32', 
                          'store_nbr': 'uint8', 
                          'family': 'category',
                          'onpromotion': 'uint16',

                          'quarter': 'uint8', 
                          'month': 'uint8', 
                          'year': 'uint16',
                          'dayofyear': 'uint16',
                          'dayofmonth': 'uint8',
                          'weekofyear': 'uint8',

                          'Friday' : 'uint8', 
                          'Monday' : 'uint8', 
                          'Saturday' : 'uint8', 
                          'Sunday' : 'uint8', 
                          'Thursday' : 'uint8', 
                          'Tuesday' : 'uint8', 
                          'Wednesday' : 'uint8', 
                          
                          'dcoilwtico' : 'float64',

                          'const': 'uint8', 
                          'trend': 'uint16',
                          })
In [ ]:
test['date'] = pd.to_datetime(test['Unnamed: 0'], format='%Y-%m-%d')
test['date'] = test.date.dt.to_period('D')
test = test.set_index('date')
test = test.drop('Unnamed: 0', axis=1)
In [ ]:
print(test.shape)
test.head(3)
(28512, 40)
Out[ ]:
id store_nbr family onpromotion Friday Monday Saturday Sunday Thursday Tuesday ... sin(6,freq=A-DEC) cos(6,freq=A-DEC) sin(7,freq=A-DEC) cos(7,freq=A-DEC) sin(8,freq=A-DEC) cos(8,freq=A-DEC) sin(9,freq=A-DEC) cos(9,freq=A-DEC) sin(10,freq=A-DEC) cos(10,freq=A-DEC)
date
2017-08-16 3000888 1 AUTOMOTIVE 0 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452
2017-08-16 3000889 1 BABY CARE 0 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452
2017-08-16 3000890 1 BEAUTY 2 0 0 0 0 0 0 ... -0.993257 -0.115935 0.796183 -0.605056 -0.154309 0.988023 -0.573772 -0.819015 0.981306 0.192452

3 rows × 40 columns

In [ ]:
stores = test['store_nbr'].unique()
families = test['family'].unique()
In [ ]:
families
Out[ ]:
['AUTOMOTIVE', 'BABY CARE', 'BEAUTY', 'BEVERAGES', 'BOOKS', ..., 'POULTRY', 'PREPARED FOODS', 'PRODUCE', 'SCHOOL AND OFFICE SUPPLIES', 'SEAFOOD']
Length: 33
Categories (33, object): ['AUTOMOTIVE', 'BABY CARE', 'BEAUTY', 'BEVERAGES', ..., 'PREPARED FOODS', 'PRODUCE', 'SCHOOL AND OFFICE SUPPLIES', 'SEAFOOD']
In [ ]:
holidays = pd.read_csv('holiday_dummies.csv')
holidays['date'] = pd.to_datetime(holidays['date'], format='%Y-%m-%d')
holidays['date'] = holidays.date.dt.to_period('D')
holidays = holidays.set_index('date')
holidays = holidays.astype('uint8')
holidays.dtypes
Out[ ]:
description_Batalla de Pichincha                      uint8
description_Black Friday                              uint8
description_Carnaval                                  uint8
description_Cyber Monday                              uint8
description_Dia de Difuntos                           uint8
                                                      ...  
description_Traslado Batalla de Pichincha             uint8
description_Traslado Independencia de Guayaquil       uint8
description_Traslado Primer Grito de Independencia    uint8
description_Traslado Primer dia del ano               uint8
description_Viernes Santo                             uint8
Length: 72, dtype: object

Funkcja wczytywanie podzbiorów treningowych¶

In [ ]:
def get_subset(store_number, family, filter,): # family jest w .query()
    current_store = pd.read_csv(f'subsets\store_nbr_{store_number}.csv', # odpowiedni format
                dtype = {'id': 'uint32', 
                        'store_nbr': 'uint8', 
                        'family': 'category',
                        'onpromotion': 'uint16',

                        'quarter': 'uint8', 
                        'month': 'uint8', 
                        'year': 'uint16',
                        'dayofyear': 'uint16',
                        'dayofmonth': 'uint8',
                        'weekofyear': 'uint8',

                        'Friday' : 'uint8', 
                        'Monday' : 'uint8', 
                        'Saturday' : 'uint8', 
                        'Sunday' : 'uint8', 
                        'Thursday' : 'uint8', 
                        'Tuesday' : 'uint8', 
                        'Wednesday' : 'uint8', 
                        
                        'dcoilwtico' : 'float32',

                        'const': 'uint8', 
                        'trend': 'uint16',
                        })
    
    current_store = current_store.query('family == @family').copy()

    current_store['date'] = pd.to_datetime(current_store['date'], format='%Y-%m-%d') # poprawka dat
    current_store['date'] = current_store.date.dt.to_period('D')
    current_store = current_store.set_index('date')

    current_store = current_store.join(holidays, how='left').fillna(0) # złączenie świąt

    # filtr sumy kumulacyjnej, z poprawką na dziedziny bez żadnej sprzedaży
    if current_store['sales'].sum() != 0:
        mask = current_store['sales'].cumsum().ne(0)
        current_store = current_store[mask]



    return current_store.query('index >= @filter') # filtr po datach z eksploracji danych
In [ ]:
for i in stores[[0]]:
    for j in families[[1]]:
        current_subset = get_subset(i, j, '2013-09-01')

current_subset
Out[ ]:
id store_nbr family sales onpromotion Friday Monday Saturday Sunday Thursday ... description_Terremoto Manabi+5 description_Terremoto Manabi+6 description_Terremoto Manabi+7 description_Terremoto Manabi+8 description_Terremoto Manabi+9 description_Traslado Batalla de Pichincha description_Traslado Independencia de Guayaquil description_Traslado Primer Grito de Independencia description_Traslado Primer dia del ano description_Viernes Santo
date
2013-09-01 433027 1 BABY CARE 0.0 0 0 0 0 1 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2013-09-02 434809 1 BABY CARE 0.0 0 0 1 0 0 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2013-09-03 436591 1 BABY CARE 0.0 0 0 0 0 0 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2013-09-04 438373 1 BABY CARE 0.0 0 0 0 0 0 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2013-09-05 440155 1 BABY CARE 0.0 0 0 0 0 0 1 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2017-08-11 2991979 1 BABY CARE 0.0 0 1 0 0 0 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
2017-08-12 2993761 1 BABY CARE 0.0 0 0 0 1 0 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2017-08-13 2995543 1 BABY CARE 0.0 0 0 0 0 1 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2017-08-14 2997325 1 BABY CARE 0.0 0 0 1 0 0 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2017-08-15 2999107 1 BABY CARE 0.0 0 0 0 0 0 0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

1445 rows × 113 columns

Funkja modeli¶

In [ ]:
from sklearn.linear_model import LinearRegression
from xgboost import XGBRegressor
In [ ]:
def hybrid_model(train, test, plot=False):
    
    ## LINEAR REGRESSION - TREND
    model_lr = LinearRegression()
    model_lr.fit(train[['trend']], train[['sales']])

    lr_fit_pred = model_lr.predict(train[['trend']])
    lr_fit_pred = pd.DataFrame(lr_fit_pred, index=train.index)


    ## ręczna kontynuacja prostej regresji liniowej, coś się krzaczyło
    step = lr_fit_pred[0][-1] - lr_fit_pred[0][-2]
    start = lr_fit_pred[0][-1] + step
    amount = 16

    manual_lr = np.linspace(start, start+(step*amount),amount,endpoint=False) 
    lr_fore_pred = manual_lr
    lr_fore_pred = pd.DataFrame(lr_fore_pred, index=test.index)

    ## detrend
    train_detrend = train.copy()
    train_detrend['sales'] =  train_detrend['sales'] - lr_fit_pred[0]


    ## XGBoostRegressor - reszta cech na zdetrendowanych danych
    model_xgb = XGBRegressor(n_estimators=25)
    X_detrend_train = train_detrend.iloc[:, 4:]
    y_detrend_train = train_detrend['sales']

    model_xgb.fit(X_detrend_train, y_detrend_train)
    detrend_test_fore = model_xgb.predict(test.iloc[:, 3:])

    test_fore = pd.DataFrame(detrend_test_fore + lr_fore_pred[0], index=test.index,).clip(0.0)
    test_fore = test_fore.rename(columns={0:'Sales'})
    test_fore['Id'] = test['id']

    ## sprawdzenie czy jest sens zaokrąglić np. sprzedaż książek to będą tylko liczby naturalne, jeśli tak - zaokrąglamy
    if train['sales'].apply(float.is_integer).all():
        test_fore["Sales"] = test_fore["Sales"].round(0)


    finished_result = test_fore.reset_index().set_index('Id').drop('date', axis=1)        

    return [lr_fit_pred, lr_fore_pred, train_detrend, test_fore, finished_result.astype('float16')]

Filtr dat¶

In [ ]:
date_filter = {
'AUTOMOTIVE' :                  '2013-01-01', # normalne 
'BABY CARE' :                   '2015-12-07', 
'BEAUTY' :                      '2013-01-01', # normalne 
'BEVERAGES' :                   '2015-05-28', 
'BOOKS' :                       '2017-01-01', # normalne 
'BREAD/BAKERY' :                '2013-01-01', # normalne 
'CELEBRATION' :                 '2015-06-01', 
'CLEANING' :                    '2013-01-01', # normalne 
'DAIRY' :                       '2013-11-07', 
'DELI' :                        '2013-01-01', # normalne 
'EGGS' :                        '2013-01-01', # normalne 
'FROZEN FOODS' :                '2013-01-01', # normalne 
'GROCERY I' :                   '2013-01-01', # normalne 
'GROCERY II' :                  '2013-01-01', # normalne 
'HARDWARE' :                    '2013-01-01', # normalne 
'HOME AND KITCHEN I' :          '2014-09-01', 
'HOME AND KITCHEN II' :         '2015-06-01', 
'HOME APPLIANCES' :             '2013-01-01', # normalne 
'HOME CARE' :                   '2015-05-04', 
'LADIESWEAR' :                  '2015-06-01', 
'LAWN AND GARDEN' :             '2016-12-03', 
'LINGERIE' :                    '2013-01-01', # normalne 
'LIQUOR,WINE,BEER' :            '2016-05-10', 
'MAGAZINES' :                   '2015-09-30', 
'MEATS' :                       '2016-10-07', 
'PERSONAL CARE' :               '2013-01-01', # normalne 
'PET SUPPLIES' :                '2015-06-01', 
'PLAYERS AND ELECTRONICS' :     '2015-06-01', 
'POULTRY' :                     '2013-11-04', 
'PREPARED FOODS' :              '2013-01-01', # normalne 
'PRODUCE' :                     '2015-06-01', 
'SCHOOL AND OFFICE SUPPLIES' :  '2013-01-01', # normalne 
'SEAFOOD' :                     '2013-01-01', # normalne 
}

Pętla prognozowania¶

In [ ]:
import random
plot_q = False # czy wykreślać i zapisywać losową część prognoz
In [ ]:
# podstawa ramki danych do dodawania kolejnych prognoz, na koniec zostaną posortowane po `ID` przed wrzuceniem wyników na Kaggle
base = pd.DataFrame(columns=['Sales', 'Id'])
base = base.set_index('Id') 
In [ ]:
for i in stores[[0]]: # na potrzeby tego zeszytu w portfolio tylko jeden sklep i jedna dziedzina zostanie zaprognozowana
    for j in families[[20]]:
        filter = date_filter[j] # która data 
        current_subset = get_subset(i, j, filter)

        current_test = test.query('store_nbr == @i & family == @j')
        current_test = current_test.join(holidays, how='left').fillna(0)

        result = hybrid_model(current_subset, current_test)

        base = pd.concat([base, result[4].astype('float16')])        

        # if random.randint(1, 100) == 1:
        # if plot_q:
        if True:
            name = j.capitalize().replace("/", " ")

            fig = go.Figure()
            fig.add_trace(go.Scatter(x=current_subset.index.to_timestamp()[-365:], y=current_subset['sales'][-365:], name=f'{j} sales')) 
            fig.add_trace(go.Scatter(x=result[3].index.to_timestamp(), y=result[3]['Sales'], mode='lines', name=f'{j} forecast'))
            fig.add_trace(go.Scatter(x=result[0].index.to_timestamp()[-365:], y=result[0][0][-365:], mode='lines', name=f'LR fit'))
            fig.add_trace(go.Scatter(x=result[1].index.to_timestamp(), y=result[1][0], mode='lines', name=f'LR fore'))

            fig.update_layout(title=f'Forecast for {name} - Store {i}', 
                            width=1200, height = 400,
                            template='plotly_dark',
                            legend=dict(
                                orientation="h",
                                yanchor="bottom",
                                y=1.02,
                                xanchor="right",
                                x=1))            
            pio.write_image(fig, f"subset images 4/{i} {name}.png")
            
            pio.write_image(fig, f"subset images cover/{i} {name}.png")



check_png(fig)
Out[ ]:
No description has been provided for this image
In [ ]:
base = base.reset_index().drop_duplicates(subset='Id').set_index('Id').sort_index()
base
Out[ ]:
Sales
Id
3000908 16.0
3002690 9.0
3004472 13.0
3006254 17.0
3008036 17.0
3009818 18.0
3011600 16.0
3013382 24.0
3015164 19.0
3016946 19.0
3018728 14.0
3020510 17.0
3022292 16.0
3024074 17.0
3025856 25.0
3027638 19.0
In [ ]:
# sprawdzenie czy rzeczywiście nie ma powtórzeń
print(base.shape)
base.index.value_counts().max()
(16, 1)
Out[ ]:
1
In [ ]:
# tu mamy ograniczoną pętlę do tylko jednego sklepu i jednej dziedziny, w pełnej petli mamy odpowiednią ilość
print(16 * 54 * 33)
print(test.shape)
28512
(28512, 40)
In [ ]:
if zapis:
    base.to_csv('my submission 4 - LR, 12 est, manual and cumsum date filter.csv')

Wnioski

Ogółem rzecz biorąc, prognozowanie wyżej przedstawioną metodą okazało się dość trafne - spośród wszystkich podejść hybryda regresji liniowej i dość prostego XGBoost'a (12 estymatorów) dała wyniki RMSLE (Root Mean Squared Logarithmic Error) równy 0.46, zajmując 88 miejsce spośród 681 uczestników. W gwoli ścisłości RMSLE liczony jest wg wzoru:

$$\sqrt{ \frac{1}{n} \sum_{i=1}^n \left(\log (1 + \hat{y}_i) - \log (1 + y_i)\right)^2}$$

gdzie:

  • $\hat y_i$ to przewidywana wartość
  • $y_i$ to rzeczywista wartość
  • $n$ to ilość obserwacji

Poniżej przedstawiono tabelę z wynikami różnych podejść:

nr podejścia Wynik Pozycja Opis
1 0.50463 122 LR + XGB, manual date filter
2 0.49233 114 LR + XGB, cumsum() date filter
3 0.50251 LR + XGB (333 estimators), cumsum() date filter
4 0.48187 108 LR + XGB (50 estimators), cumsum() date filter
5 0.47315 98 LR + XGB (25 estimators), cumsum() date filter
6 0.46249 88 LR + XGB (12 estimators), manual and cumsum() date filter
6 0.46807 LR + XGB (5 estimators), manual and cumsum() date filter
7 0.47727 Ridge $\alpha=1$ + XGB (25 estimators), manual and cumsum() date filter
8 0.46909 Ridge $\alpha=3$ + XGB (12 estimators), manual and cumsum() date filter
9 0.46909 Ridge $\alpha=9$ + XGB (12 estimators), manual and cumsum() date filter

Okazało się, że zaawansowany model w hybrydzie powinien być dość prosty - najlepsze wyniki przy 12 drzewach (estymatorach) w XGBoost. Regresja grzbietowa (Ridge Regression) zmniejszająca dopasowanie do danych treningowych nie przyniosła lepszych efektów.

Zalecenia jak usprawnić model:

  • dokładniejsza eksploracja danych - to jest niestety trudne w przypadku podejścia, gdzie każdy sklep i jego dziedziny sprzedaży modeluje się osobno
  • lepsze przyjrzenie się świętom - prawdopodobnie wrzucenie wszystkich świąt niepotrzebnie powoduje szumy, na podstawie eksploracji można by zaznaczyć kilka świąt / wydarzeń - Wielkanoc, początek roku szkolnego (sprzedaż artykułów biurowych i szkolnych), Święta Bożego Narodzenia, Nowy Rok, istotniejsze święta narodowe Ekwadoru
  • zastosowanie innych modeli np. regresji LASSO, Elastic Net, sieci neuronowych
  • zastosowanie zupełnie innego podejścia - np. modelowania ARIMA specjalnie przeznaczonego do serii czasowych
  • zastosowanie rekursywnej metody predykcji:
    • każda pojedyncza predykcja (sprzedaż danego dnia, danej dziedziny w danym sklepie) zostaje dołączona do zestawu treningowego i kolejno następuje ponowne nauczanie modelu i kolejna pojedyncza predykcja
    • zaletą tej metody mogą być celniejsze wyniki, z drugiej strony wczesny błąd w predykcjach będzie "niesiony" z każdą kolejną predykcją,
    • metodę tę można poprawić dodając lagi wartości, którą chcemy przewidzieć (tu już pojawiają się zagadnienia z ekonometrii), lagi to wartości z przeszłości
    • główną wadą metody rekursywnej są wysokie wymagania obliczeniowe